home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / nrpas13.zip / MPROVE.DEM < prev    next >
Text File  |  1991-04-29  |  2KB  |  64 lines

  1. PROGRAM d2r5(input,output);
  2. (* driver for routine MPROVE *)
  3. CONST
  4.    n=5;
  5.    np=n;
  6. TYPE
  7.    glnarray = ARRAY [1..n] OF real;
  8.    glindx = ARRAY [1..n] OF integer;
  9.    glnpbynp = ARRAY [1..np,1..np] OF real;
  10. VAR
  11.    glinext,glinextp : integer;
  12.    glma : ARRAY [1..55] OF real;
  13.    d : real;
  14.    i,idum,j : integer;
  15.    a,aa : glnpbynp;
  16.    b,x : glnarray;
  17.    indx : glindx;
  18.  
  19. (*$I MODFILE.PAS *)
  20. (*$I RAN3.PAS *)
  21.  
  22. (*$I LUDCMP.PAS *)
  23.  
  24. (*$I LUBKSB.PAS *)
  25.  
  26. (*$I MPROVE.PAS *)
  27.  
  28. BEGIN
  29.    a[1,1] := 1.0; a[1,2] := 2.0; a[1,3] := 3.0; a[1,4] := 4.0; 
  30.    a[1,5] := 5.0; a[2,1] := 2.0; a[2,2] := 3.0; a[2,3] := 4.0;
  31.    a[2,4] := 5.0; a[2,5] := 1.0; a[3,1] := 1.0; a[3,2] := 1.0;
  32.    a[3,3] := 1.0; a[3,4] := 1.0; a[3,5] := 1.0; a[4,1] := 4.0;
  33.    a[4,2] := 5.0; a[4,3] := 1.0; a[4,4] := 2.0; a[4,5] := 3.0;
  34.    a[5,1] := 5.0; a[5,2] := 1.0; a[5,3] := 2.0; a[5,4] := 3.0;
  35.    a[5,5] := 4.0;
  36.    b[1] := 1.0; b[2] := 1.0; b[3] := 1.0; b[4] := 1.0; b[5] := 1.0;
  37.    FOR i := 1 to n DO BEGIN
  38.       x[i] := b[i];
  39.       FOR j := 1 to n DO BEGIN
  40.          aa[i,j] := a[i,j]
  41.       END
  42.    END;
  43.    ludcmp(aa,n,np,indx,d);
  44.    lubksb(aa,n,np,indx,x);
  45.    writeln;
  46.    writeln('Solution vector for the equations:');
  47.    FOR i := 1 to n DO write(x[i]:12:6);
  48.    writeln;
  49. (* now phoney up x and let mprove fit it *)
  50.    idum := -13;
  51.    FOR i := 1 to n DO BEGIN
  52.       x[i] := x[i]*(1.0+0.2*ran3(idum))
  53.    END;
  54.    writeln;
  55.    writeln('Solution vector with noise added:');
  56.    FOR i := 1 to n DO write(x[i]:12:6);
  57.    writeln;
  58.    mprove(a,aa,n,np,indx,b,x);
  59.    writeln;
  60.    writeln('Solution vector recovered by mprove:');
  61.    FOR i := 1 to n DO write(x[i]:12:6);
  62.    writeln
  63. END.
  64.